ggplot2绘制曼哈顿图示例
备注：在这篇文献中，使用R包edgeR中的方法鉴别显著富集的OTUs。材料方法部分截图如下。
ggplot2绘制微生物组分析中的曼哈顿图示例
作图数据说明
ggplot2作图示例
首先加载R包并读取作图数据。
library(ggplot2)#读取数据otu_stat <- read.delim('otu_sign.txt', sep = '\t'原始表格里面OTUs的顺序根据id排列的，而我们看到上述文献曼哈顿图中的OTUs则是根据其所属分类排列的，因此我们需要首先根据OTUs所属分类排个序（文献中按目分类，示例数据为了方便演示直接按门分类了），便于在后续作图时按分类水平分别呈现OTUs。
这里排序的时候直接根据门分类名称的首字母排序了，按门分类排序后，对于每个微生物门中的OTUs顺序，则使用默认顺序。如果你比较讲究，想按某种顺序排列，自定义调整即可。
#门水平排序，这里直接按首字母排序了otu_stat <- otu_stat[order(otu_stat$phylum), ]otu_stat$otu_sort <- 1:nrow(otu_stat) #其它自定义排序，例如你想根据 OTU 数量降序排序#phylum_num <- phylum_num[order(phylum_num, decreasing = TRUE)]#otu_stat$phylum <- factor(otu_stat$phylum, levels = names(phylum_num))#otu_stat <- otu_stat[order(otu_stat$phylum), ]#otu_stat$otu_sort <- 1:nrow(otu_stat)然后作图就可以了。为了更好理解作图过程，暂不考虑调整背景细节，先将散点绘制出，得到曼哈顿图的初始轮廓。
#初始样式p <- ggplot(otu_stat, aes(otu_sort, -log(p_value, 10))) + #geom_point(aes(size = abundance, color = phylum, shape = enrich)) + #点的大小表示 OTUs 丰度，颜色表示分类，形状表示是否为富集 OTUsscale_size(range = c(1, 5))+ #点大小范围scale_shape_manual(limits = c('sign', 'no-sign'), values = c(16, 1)) + #点形状，实心点为富集 OTUs，空心点为非富集 OTUs labs(x = NULL, y = '-log10(P)', size = 'relative abundance (%)', shape = 'enriched') + #坐标轴标题、图例标题等theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + #主题调整，去除背景线等guides(color = 'none') + #隐藏颜色图例（OTUs 分类）geom_hline(yintercept = -log10(0.01), color = 'gray', linetype = 2, size = 1) #在 p=0.01 处的 -log10 位置（即 y = 2）绘制条横虚线，表示 p<0.01 为差异 OTUs 的判定标准之一 p初始的结果如下，大致上曼哈顿图的主体部分就出来了（先不管x轴刻度，放在下一步中说）。这里我没有再去设置颜色分类，使用的ggplot2的默认颜色，大家有需要可以自行调整。
然后，我们在横坐标的对应位置处，将微生物门类群的标签添加上。
由于横坐标是根据（按微生物门分类）排序后的OTUs绘制的，因此我们首先需要计算各门类群所包含OTUs的排序序号的中值位置，即得到标签坐标。在计算中值位置，用于得到各微生物门标签坐标的同时，也计算各门类群所包含OTUs的排序序号的边界位置，在后续将根据边界坐标绘制矩形区域，使区域区分更明显。
#计算 x 轴标签、矩形区块对应的 x 轴位置phylum_num <- summary(otu_stat$phylum) phylum_range <- c(0, phylum_num[1])phylum_name <- phylum_num[1] / 2for (i in 2:length(phylum_num)) { phylum_range[i+1] <- phylum_range[i] + phylum_num[i] phylum_name[i] <- phylum_range[i] + phylum_num[i] / 2}然后在计算好的横坐标位置处，添加微生物门分类刻度标签。
#添加 x 轴刻度标签p <- p +scale_x_continuous(breaks = phylum_name, labels = names(phylum_num), expand = c(0, 0)) +theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) p然后再添加矩形区域。这里使用ggplot2中的annotate()函数绘制额外的矩形，annotate()可用于在ggplot2中的指定坐标位置处添加额外的形状、箭头、文字标签等，详情使用?annotate查看帮助。方便起见，矩形的颜色就直接交替绘制为不同深度的灰色了。
#添加矩形区块，交替绘制为不同深度的灰色#由于先绘制的散点，后绘制的矩形，即矩形图层在上，故需要设置很高的透明度才可for (i in 1:(length(phylum_range) - 1)) p <- p + annotate('rect', xmin = phylum_range[i], xmax = phylum_range[i+1], ymin = -Inf, ymax = Inf, alpha = 0.1, fill = ifelse(i %% 2 == 0, 'gray60', 'gray40')) p这样一个与文献中风格类似的曼哈顿图就算搞定了。
上文为了更好地演示过程，就先绘制了散点后，再绘制了矩形区域，这样矩形区域的图层就位于了散点的上方。尽管我们对矩形区域设置了透明度，尽可能降低了矩形区块对散点的覆盖，但是相比之下，肯定是先把矩形绘制好，再添加散点的效果会更佳。
##推荐先调整好背景（矩形区块），再绘制散点，这样散点位于矩形区块图层的上面，更清晰#背景、标签布局p <- ggplot(otu_stat, aes(otu_sort, -log(p_value, 10))) +labs(x = NULL, y = '-log10(P)', size = 'relative abundance (%)', shape = 'significantly enriched') +theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) +scale_x_continuous(breaks = phylum_name, labels = names(phylum_num), expand = c(0, 0)) +theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) #这样也没必要再为矩形设置透明度了for (i in 1:(length(phylum_range) - 1)) p <- p + annotate('rect', xmin = phylum_range[i], xmax = phylum_range[i+1], ymin = -Inf, ymax = Inf, fill = ifelse(i %% 2 == 0, 'gray95', 'gray85')) #最后绘制散点p <- p + geom_point(aes(size = abundance, color = phylum, shape = enrich)) +scale_size(range = c(1, 5))+scale_shape_manual(limits = c('sign', 'no-sign'), values = c(16, 1)) + guides(color = 'none') +geom_hline(yintercept = -log10(0.01), color = 'gray', linetype = 2, size = 1) p #输出图片至本地ggsave('manhattan.pdf', p, width = 10, height = 5)ggsave('manhattan.png', p, width = 10, height = 5)上述各步命令详解可参见上文，除了顺序，基本没做细节的更改。
输出图片结果如下，这样我们可以选择使用更深色的矩形颜色，区域区分更明显。
后期修图调整
对于输出在本地的结果，可以考虑使用AI、PS等工具后期调整，达到更好的效果。如下示例，我在AI（Adobe Illustrator）中简单地调整了输出的pdf文件（矢量图格式，比png位图更方便调整）中微生物门分类标签的位置，图例圆圈的位置、大小等，这样我们的作图结果就和文献中的样式更接近了。其实很多细节部分的调整，后期修图比在R中调试参数更为节省时间，如上文中的矩形绘制，个人就感觉其实在AI中添加比在R中绘制方便。
方法大致就这样吧。当然，模仿文献中图形的样式不是重点，本文的目的还是参照文献对ggplot2作图方法作个简介，更多更美观的可视化样式，还得细致去琢磨了。
ggplot2绘制GWAS分析中的曼哈顿图示例
这部分就不详细说了，方法基本同上，就尽管数据类型换了，但实际上数据的排列方式没啥太大改变。况且，关于GWAS曼哈顿图的绘制方法，网上的相关示例很多，所以咱也没必要多此一举了……因此只简单演示下。
示例数据就直接使用qqman包中的gwasResults数据集了。
#直接使用了 qqman 包中的 gwasResults 数据集gwasResults <- qqman::gwasResultsfix(gwasResults) #查看这个数据集还是很好理解的，为人类的22条常染色体中，与某性状相关的SNP位点信息了。SNP，每一个SNP位点的序号；CHR，该SNP位点位于第几条染色体上；BP，该SNP位点在该染色体中的序号；P，关联程度的显著性p值。
在qqman包中，使用manhattan()等函数即可直接绘制曼哈顿图，更换其他数据同样可行。详见?manhattan帮助说明。
#qqman 包 manhattan()绘制曼哈顿图library(qqman)manhattan(gwasResults, col = c('gray90', 'gray80')) #使用两种灰度的颜色交替绘制各染色体中的 SNP 位点但此处我们不再过多介绍它们（想了解的话，网上相关教程挺多的），而是使用ggplot2绘制曼哈顿图。一个示例如下，作图参数大致同上，不再详说。颜色设置也直接使用ggplot2默认配色。
#这里借助 doBy 中的 summaryBy() 计算各染色体中 SNP 序号的中值位置，以放置染色体标签library(doBy) gwasResults$SNP1 <- seq(1, nrow(gwasResults), 1)gwasResults$CHR <- factor(gwasResults$CHR, levels = unique(gwasResults$CHR))chr <- summaryBy(SNP1~CHR, gwasResults, FUN = median) #ggplot2 作图p <- ggplot(gwasResults, aes(SNP1, -log(P, 10), color = CHR)) +theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'), panel.background = element_rect(fill = 'transparent')) +geom_point(aes(color = CHR), show.legend = FALSE) +labs(x = 'Chromosome', y = expression(''~-log[10]~'(P)')) +scale_x_continuous(breaks = chr$SNP1.median, labels = chr$CHR, expand = c(0, 0)) +geom_hline(yintercept = c(5, 7), color = c('blue', 'red'), size = 0.5) p #保存ggsave('GWAS.pdf', p, width = 10, height = 4)ggsave('GWAS.png', p, width = 10, height = 4)